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The ommatidium of the butterfly's afocal apposition eye exhibits angular performance that can only be 
achieved by transforming the diffraction pattern of its corneal lens into the fundamental mode of its 
rhabdom waveguide. A graded index model of the ommatidium has been proposed and verified but the 
efforts to extract the transformation's underlying physics from it have been hindered by its extreme 
complexity. Here we numerically investigate the ommatidium model and reveal that the current model, 
involving only the graded index distribution, is insufficient for the transformation. We further find that 
adding spatially varying absorption to the existing model dramatically improves its transformation 
performance, producing near-perfect mode matching with overlap integral exceeding 0.96. Such a 
combined action of spatially varying index and absorption for microscale mode transformation is new to 
researchers in optics and biology and will benefit both disciplines. 

I nsects' compound eyes achieve their wide field-of-view using hemispherically arrayed optics 1 . In the apposi- 
I tion type, one corneal lens (CL) and one rhabdom (Rh) form a lens- waveguide pair which, along with other 
I parts between the two, is collectively called an ommatidium (Fig. la). For optical isolation between ommatidia, 
the level of light coupled into Rh depends critically on the incidence angle 0i. The angular response typically peaks 
at normal incidence (6i = 0) and tapers off as 0i increases. The acceptance angle 0 A marks the full-width half- 
maximum point of the response curve 1 " 3 . In the superposition type, in contrast, each CL collects rays over a wide 
angular range. Then the crystalline cone (CC) collimates them for redistribution over multiple photoreceptors 
(Fig. lb) 1 ' 4 . 

Nilsson et at. found a third type in butterflies 5 . Structurally, it is an apposition eye comprising a large number of 
ommatidia. But it also resembles the superposition eye in collimating its output using the cone stalk (CS) as the 
collimation lens (Fig. lc). These mixed characteristics positioned this third type, called afocal apposition eye for 
its perfect collimation, as the link between the apposition and superposition eyes 5 ' 6 . 

At question was how an apposition eye came to employ collimation which does not give clear advantages while 
requiring very complex CS architectures. For instance, with its unique index distribution, Heteronympha meropes 
CS produces 0.2 mega-diopters of magnification from just —10 urn in length, making itself "the most powerful 
known to man. 5 " 

Nilsson et al. provided a clue by revealing another role of CS 7 . As mentioned above, ommatidial isolation 
requires reductions in 0 A . Optically, 0 A minimizes when the electric field pattern hitting the Rh waveguide 
entrance exactly matches that of the waveguide's fundamental mode (LP 01 ). However, CL generates Airy- like 
patterns that are different from LP 01 . In typical apposition eyes, the mismatch broadens 0 A by —15% from its 
theoretical minimum 6 ' 7 . In contrast, the afocal apposition eye's 0 A turns out to match the theoretical minimum, 
indicating that CS transforms the Airy-like pattern into LP 01 . Collimation is just a part of the mode transforma- 
tion process. 

The search for the transformation mechanism has produced a graded index model of the ommatidium 8 and 
multiple hypotheses. One interpreted CS as a waveguide mode coupler 7 . Since Airy pattern closely resembles the 
sum of the two lowest axis -symmetric modes, LP 01 and LP 02 , CL practically excites both modes into CS. Totally 
transferring LP 0 2's mode power towards LP 01 will realize the LP 01 transformation. Another suggested that CS 
functions as a mode mixer which synthesizes LP 01 from the local modes of the tapered, index-graded CS 8 . So far, 
theoretical efforts to verify the hypotheses have been thwarted by the model's extreme complexity. 

Here we analyze the LP 0 i transformation mechanism numerically. We find that the existing graded index 
model is insufficient for the task. Complementing it with spatially varying absorption, however, has directly led to 
successful LP 01 transformations. Such a combined action of non-uniformly distributed index and absorption 
gives the afocal ommatidium a new interpretation as an absorption-assisted mode shaper, rather than a coupler or 
mixer. The new findings will enrich the knowledgebase of insect optics 1 and bio-inspired optics 910 alike. 
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Figure 1 | Features of afocal apposition eye and its model. (a,b) Schematic diagrams of (a) apposition and (b) superposition compound eyes. The former 
consists of CL-Rh pairs. Each pair maximally accepts only the normally impinging ray with the level of acceptance tapering off as the incidence angle 
increases. The full-width half-maximum point of the curve marks the acceptance angle 6 A . In the latter, in contrast, Ph collects rays from multiple CLs. 
Light collimation by CC is essential for propagation across CZ. (c) The ommatidium of the afocal apposition eye. CL focuses light inside CC, not at the Rh 
entrance as in (a). Similar to the superposition eye, its CS collimates the defocusing rays. (CL: corneal lens, CC: crystalline cone, CS: cone stalk, Rh: 
rhabdom, Ph: photoreceptors, CZ: clear zone). 



Results 

Testing mode coupling/mixing hypotheses. With the numerical 
simulation code and ommatidium model set up as described in 
Methods, we first simulated the propagation of a plane wave 
through the afocal ommatidium to test the existing hypotheses. 
Figures 2a and b show that the simulation correctly reproduced the 
experimentally observed mid-CC focal spot formation (z ~ 58 |im). 
Despite the non-uniform index distribution in CC, the focal field 
pattern closely matched Airy pattern, both in amplitude and phase 
(Fig. 2c). Existing studies hypothesized that this Airy-like pattern 
would evolve into LP 01 through mode coupling or mixing as it 
approaches the Rh entrance at z = 68 |im 7 ' 8 . 

However, the amplitude and phase profiles obtained at that point, 
shown in Fig. 2d, exhibited high-level mismatches from those of 
LP 0 i. The side-lobes of the Airy- like pattern persisted throughout 
CS, leaving severe wiggles to the final amplitude profile. The step-like 
phase distribution in Fig. 2c has been quasi-linearized as shown in 
Fig. 2d but it is still insufficient to match the flat wavefront of LP 01 . 
Only the width of the field pattern matched that of LP 01 reasonably at 
the Rh entrance. 

To quantify the level of mode matching, we introduced three 
metrics in Methods: the shape error 5 S , phase error 5 p , and the 
overlap integral rj which considers both shape and phase. The inset 
in Fig. 2b shows the computed rj as a function of z. Even within CS, r\ 
increases very slowly, only to reach 0.77 at its end. The values of 5 S 
and 5 p are also large, reaching 0.0603 and 0.6068 at the Rh entrance, 
respectively. Simulations performed at different wavelengths with 
slightly varied structural parameters, including the CP thickness, 
CC length, and CS index profile, all produced similar mismatch. It 
indicates that mechanisms relying solely on the graded index distri- 
bution of CS, such as mode coupling or mixing, are insufficient to 
transform the Airy-like pattern into LP 0 i- 



Impact of distributed apodization. The intensity distribution in 
Fig. 2a shows that the wiggles in the final amplitude profile 
evolved mainly from the outermost fringe of the CL-focused light. 
So we attempted to suppress the wiggles by blocking the fringe with 
an annular absorptive region. It is analogous to the apodization 
process in Fourier optics which removes the peripheral, high 
spatial frequency portion of a diffraction pattern 11 . Removal of the 
wiggles alone, however, does not guarantee LP 0 i transformation. It 
hinges on whether the subsequent geometry and index distribution, 
already given by Nilsson et at. 5,6 and Pask et al. 8 , can correctly modify 
the amplitude and phase of the apodized wave. 

We explored its feasibility with iterative simulations. We assumed 
that the absorption is mild yet broadly distributed and implemented 
it by making the portion of CC with p > a and z ^ z 0 weakly 
absorbing in our numerical model, as shown in Fig. 3a. We empir- 
ically set the electric field attenuation at 1% per 50 nm or 
0.873 dB-|im _1 . Given the existence of light-absorbing pigments 
around the ommatidium, this type of distributed absorption in CC 
seems plausible 1213 . Figure 3b shows the intensity distribution from 
an exemplary simulation. The absorptive region apparently blocked 
most of the outermost fringe. 

The apodization impacted LP 01 transformation immediately. 
Figure 3c shows the amplitude and phase profiles obtained at 
z = 67 A |im, right in front of the Rh entrance, with (z G , a) = (33, 
2.4) [im. Comparison of the amplitude profile with that of LP 0 i 
reveals that the apodization was very effective not only in removing 
the wiggles but also in reshaping the overall mode profile into that of 
the single-peaked LP Ql . The computed shape error 5 S was 0.0216, a 
factor of 3 reduction from that of the unapodized one. 

It is important to point out that the distributed absorption greatly 
improved the phase profile matching as well, even though it cannot 
directly affect phase. At the same z point (67.4 |im), the phase error 
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Figure 2 | Testing mode coupling/mixing hypotheses, (a) and (b) show the intensity distribution and the change in its center axis intensity resulting from 
the simulated propagation of a plane wave through the ommatidium model in Methods. Both confirms the mid- CC focusing by CL near z = 58 (am. The 
inset in (b) shows the change in the overlap integral rj between z = 55 and 68 urn. (c) The amplitude and phase profiles at the focal spot (z = 58 (am). They 
agree well with those of Airy pattern except for the peripheral, low-intensity part (p > 5 (am), (d) The amplitude and phase profiles obtained at 
the Rh entrance (z = 68 (am), however, exhibit significant mismatches from the Rh waveguide's LP 01 mode, indicating that the hypothesized mode 
coupling or mixing mechanism is insufficient for transforming an Airy- like pattern into LP 01 . 



5 p was 0.0178, a factor of 34 improvement from the unapodized one. 
The inset of Fig. 3c visualizes the flatness of the wavefront across the 
radial extent of CS (p < 1.3 (am). We attribute this improvement 
again to the apodization which suppressed the peripheral portion 
with its phase quadratically retarded by CL. 

The results of iterative simulations shown in Fig. 4 further reveal 
that such large 5 p improvements near the Rh entrance is possible only 
within a narrow window with 33 < z Q < 36 (am and 2.1 < a < 
2.5 (am, reaffirming the need for a precise and complete elimination 
of the outermost fringe. A simultaneous reduction of 5 S within the 
window, however, was hindered by the persistent bulging near p = 
2 (am (Fig. 3c). Since the level of shape matching within the radial 
extent of CS remains excellent, the bulging must be due to the 



absence of absorption outside the terminal portion of CS, which is 
evident in Fig. 3a. 

Refinement by conformally extended absorption. So we extended 
the absorption region as shown in Fig. 5a. Starting from z = z t , the 
exterior of CS was made absorptive all the way to the original 
apodization region. A thin buffer layer with thickness t was 
installed for fine control. For simplicity, we set the electric field 
attenuation within this conformally extended absorption region 
identical to that of the apodization region, i.e., 0.873 dB-(am _1 . 
Then we performed iterative simulations involving all four 
parameters: z Q (33 ~ 36 (am), a (2.1 ~ 2.5 |am), z t (65 ~ 67 (am), 
and t (0 ~ 0.5 (am) and again observed additional improvements in 
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Figure 3 | Absorptive apodization and its impact, (a) The setting of the hypothesized absorptive apodization. Low-level absorption (0.873 dB (am" 1 ) was 
implemented in the region with z > z Q and p > a. (b) Intensity distribution obtained from an apodized propagation simulation with (z 0 , a) = (36, 
2.3) (am. Blocking of the outermost fringe is evident, (c) The amplitude and phase profiles obtained at z = 67.4 (am at which the phase error 5 p becomes 
minimized. The dotted red curves represent LP 0 i which also appeared in Fig. 2 d. The bulging in the amplitude profile for p ~ 2 (am persisted for other 
values of (z G , a), necessitating additional mode shaping. 



LP 01 transformation. We present the results obtained with (z G , a) = 
(36, 2.4) (am as examples. 

Figure 5b shows the amplitude profile obtained at z = 67.9 (am 
with the value of t fixed at 0.35 (am. When z t = 65 (am, the bulging 
near p = 2 (am disappeared completely, enabling the amplitude 
profile to overlap that of LP 01 almost perfectly. The corresponding 
5 S = 0.0049, a factor of 5.3 reduction from the Fig. 3c result. The 
absolute difference between the amplitude profiles of the simulation 
results and LP 01 , plotted in the inset, shows that the match improves 
increasingly as the conformally extended absorption region becomes 
longer. Figure 5c shows the phase profile at the same z-point as a 
function of z t . Just as the apodization did, this conformally extended 
absorption further reduced the phase error values by flattening the 
wavefront. When z t = 65 (am, it became near-flat within the radial 
extent of CS. The simultaneous minimization of 5 p and 5 S enabled the 
overlap integral y\ to reach 0.9646 at the Rh entrance. 

The minimum phase and shape error values (5 p5min and 5 Sjmin ) 
plotted in Fig. 5d show that simultaneous 5 p - 5 S minimization is 
limited to a narrow window of (z t , t) combinations. The values of 
the optimal (z t , t) combinations were also different for different (z oy 



a) combinations. So, we extracted from our iterative simulation 
data the combinations of (z Q , a, z t , t) which minimize both 5 p 
and 5 S at z > 67.5 (am, i.e., within 500 nm of the Rh entrance, 
while maintaining 5 S < 0.01. Table 1 lists them along with their 
performance metrics. It shows that the optimal values for z 0 fall 
within a 3 (am range from the CC starting point, with the majority 
at z Q = 35 (am. The optimal values for the apodization pupil radius 
a are unanimously at or near 2.5 (am, which is approximately 42% 
of the CC radius at its starting point. The optimal z t and t values 
exhibit only weak variations within 1 — 3 (am and 0.2 ~ 0.4 (am 
range, respectively. As shown in Fig. SI in Supplementary 
Information, these configurations are effective over a wide spectral 
range spanning 400 to 600 nm. 

Since the proposed LP 01 transformation relies on absorption, it 
inevitably reduces the mode power. For optimal configurations in 
Table 1, they can be 73 ~ 77% of the initial value. As shown in Fig. 5e, 
the initial apodization claims ~5 dB and the extended absorption 
another ~1 dB. The overall loss of 6 dB may seem high in power 
budget point of view. Removing the collimating CS and the absorpt- 
ive mode transformation altogether, however, results in a situation in 
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Figure 4 | Parametric search for simultaneous minimization of 5 S and 5 p . 

(a) Solid curves: z-positions at which the phase and shape errors minimize 
for the given (z Q , a). Lack of their intersection indicates that the two do not 
minimize simultaneously. Dashed curves: z min)er ror values obtained by 
limiting the range for the shape error evaluation within p < 1.3 (am, which 
enabled simultaneous minimization of 5 S and 8 p . (b) 5 P)m i n values as a 
function of (z OJ a). In general, 5 p grows rapidly as z Q moves away from CP. 
(c) Solid curves: 5 Sjm i n values obtained when z Q = 33 ~ 36 (am, i.e., within 
3 (am from the ending point of CP. Overall, 5 S values grow monotonically 
with a. Dashed curves: 5 S assessed only within p < 1.3 (am. Now they 
exhibit definite dips, suggesting that the shape error accumulates mostly in 
the region with p > 1.3 (am, i.e., the area not covered by the distributed 
absorption shown in Fig. 3a. 

which the exact Airy pattern tries to excite LP 0 i- The corresponding 
coupling efficiency is 14 : 

where x = 7.34-(r 0 /r A ) 2 , r A is the Airy spot radius defined as 
O.61-/cl"^o/( w cl"0cl) an d r 0 is the modal field spot size for 
Gaussian-approximated LP 01 which is —1.7 (am for the present Rh. 



The resulting rj coup i e is —6.7 dB. Table 1 shows that the present 
model incurs slightly lower loss while adding the advantage of LP 01 
transformation and, hence, 0 A minimization. This justifies the use of 
absorption for mode transformation. The last column of Table 1 also 
indicates that initiating the apodization early, i.e., near the CP/CC 
interface, can increase rj slightly but also incurs significant penalty to 
the mode power. 

Discussion 

In conclusion, we utilized numerical simulations to test the existing 
hypotheses on how the afocal ommatidium transforms the CL-gen- 
erated Airy-like pattern into LP 01 and push its acceptance angle 
towards the theoretical minimum. Our initial results showed that 
the hypothesized mechanisms relying solely on graded index distri- 
bution, such as mode coupling or mixing, cannot accomplish the 
LP 01 transformation. We then found that complementing the exist- 
ing graded index model with spatially varying absorption can dra- 
matically increase the level of mode matching. Interestingly, the 
added absorption improved not only the shape matching, which is 
under direct influence of loss, but also the wavefront flatness. This 
wavefront-flattening operation seems to develop into the collimation 
step of the superposition eye. Such a combined action of distributed 
index and absorption for microscale mode transformation is novel 
with no previous report. Iteratively, we identified a window in para- 
meter space capable of realizing LP 0 i transformation with overlap 
integral exceeding 96%. The 6 dB power penalty arising from the use 
of absorption can easily be justifiable given that the purely refractive 
model also imposes a comparable level of loss due to the Airy-to-LP 01 
mode mismatch even though it cannot bring in the added advantage 
of 0 A minimization predicted for the absorption-assisted model. 

It needs to be emphasized that the LP 0 i transformation in this 
work was achieved with minimal modifications made to the already 
verified model of the afocal ommatidium. In biological ommatidia, 
those modifications, i.e., distributed absorption within and around 
CC, can be realized with light absorbing pigments infiltrating and 
surrounding the ommatidia, respectively. In soft biological struc- 
tures, achieving and maintaining the (am-scale dimensional precision 
for optimal LP 0 i transformation seems challenging. Given the exist- 
ence of pigments capable of migrating around the ommatidium for 
dark/light adaptations, we can envision a possible adaptive control 
scheme which adjusts the pigment distribution to maximize the 
mode matching. 

These findings will interest both optical engineers and biologists. 
The former may get excited at the prospect of realizing highly min- 
iaturized mode converters that are coupled with microscale focusing 
lenses and/or formed in-fiber. The recent trend in optical commun- 
ication to exploit various modes supported by multiple mode fibers 21 
indeed necessitates new paradigms in optical mode conversion 22,23 
and the new absorption-assisted scheme will certainly contribute to 
the efforts. The possible involvement of absorption in the operation 
of insects' ommatidia would motivate the latter to examine addi- 
tional roles of absorptive pigments. Numerically reproducing the 
angular response of the afocal apposition ommatidium using a fully 
3 -dimensional model will give the present work an ultimate valid- 
ation. Here we emphasize that the optimal configurations found in 
this work are by no means ultimate or exhaustive. They can be 
further refined once the currently unknown refractive index and loss 
distributions around CC become available and the possibility of 
finding purely refractive models capable of performing the required 
transformation, as suggested in Ref. 23, should not be ruled out 
either. 

Methods 

Numerical simulation of wave propagation. To efficiently simulate the optical wave 
propagation in the ommatidium, we utilized the finite-difference beam propagation 
method (FD-BPM) 15 . Specifically, we chose a version of FD-BPM customized for 
scalar wave propagation in an axially symmetric cylindrical space 16 " 20 and 
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Figure 5 | Impact of conformally extended absorption, (a) The setting of the conformally extended absorption which starts from z = z t with a f um-thick 
buffer. (b,c) The amplitude and phase profiles obtained at the Rh entrance when t = 0.35 um. LP 01 profiles are also superimposed for comparison, (b) For 
= 65 (am, i.e., 3 um-long conformal absorption, the output becomes indistinguishable from LP 0 i- The inset shows how the increase in the conformal 
absorption reduces the mismatch in the output amplitude profile, especially near p = 2 (am. (c) The conformal absorption also improves the wavefront 
flatness within the radial extent of CS. (d) The minimum phase and shape errors for various (z t , t) combinations. Note that 5 p and 5 S minimize 
simultaneously near (zt, t) = (65 ~ 66, 0.2 ~ 0.35) (am. (e) The change in the mode power along the propagation through CC and CS. CP/CC: The 
interface between CP and CC. 



implemented it with Matlab. Azimuthally changing modes, such as LP X i or LP 2 i, were 
inherently excluded. More universal methods exist for vectorial, azimuthally varying 
wave propagation but their computational burden was prohibitive for iteratively 
simulating a structure measuring over 15 um in radius and 70 um in length. 
Throughout our work, the radius of the simulation space was 20 um. We used the 
transparent boundary condition to suppress spurious reflections from the edges of 
computation domain. The spatial resolution was fixed at 50 nm in both radial and 



axial directions. A full propagation typically took 15 — 30 seconds on computers with 
a xeon processor and 4 GB RAM. We validated the FD-BPM code by simulating the 
formation of Airy pattern by a spherical lens (Supplementary Figure S2). The 
wavelength X Q was varied between 400 and 600 nm. 

Heteronympha merope ommatidium model. Most structural and optical 
parameters in the ommatidium model were from Nilsson et al. 5 ' 6 and Hateran et al 7 : 
CL radius: acL = 12.25 um, CL thickness: 23 um, CL radius of curvature: 20 um, 
n C L = 1-52./cl ~ 54 um. CP was modeled as a 10 um-thick plate with n C p = 1.34. 
The Rh parameters are also close to those of Nilsson et al. 5 ' 6 : Rh diameter = 2.07 ~ 
3 um, «Rh_amb = 1-35, and = 1.36 — 1.38. For the graded index profile and 
dimensions of CC and CS, we adopted the model proposed by Pask et al. 8 : n cc 2 = 
n Q 2 -[l + 2-A-(G! 2 - G 2 )/P 2 ] where n Q = 1.4, A = 131.5, 0! = 0.1 rad. R and 0 were 
computed at each position within the CC. With the length of CC set to 35 um 8 , its 
radius decreases linearly from 6.1 to 1.3 um as it approaches Rh. Pask et al. verified 
this model 8 by reproducing the experimental optical magnification data of Nilsson et 
al. 5 Supplementary Figure S3 shows the numerical model schematically. 

Quantification of mode matching. To quantify how closely a simulated electric field 
distribution matches that of LP 0 i, we utilized three metrics. The shape error 8 S is 
defined as: 



Os = y|(|^| 2 -|P, e/ | 2 ) 2 ^ /j|E re/ | 2 <Js (2) 

where E sim and E re f represent the electric field distributions obtained from the 
simulation and the exact LP 0 i mode solution, respectively. Since our FD-BPM 
presumes total axis -symmetry, the cross-sectional integration was actually performed 
only along the radial direction. Two different integration ranges were used. For 



Table 1 | Configurations for LPq] transformation and results. (z Q , 
a, z t , fj combinations producing simultaneous 5 p - 5 S minimization 
at z > 67.5 um and 5 S < 0.01 . The corresponding 5 S and 5 p 
values can be found in the second column. The ratio between the 
final and initial mode power and the overlap integral rj are given 
in the third column 
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overall shape error, the intensity distributions within p = 0 ~ 7 urn range was 
compared. When comparing only the center portion of the intensity distribution, the 
range was reduced to p = 0 ~ 1.3 um. Since LP 01 exhibit a flat phase distribution 
across its mode profile, we assessed the level of match in the phase with the phase 
error § p which is the standard deviation of the simulated phase response. Because CS 
manipulates only the center portion of the wave, 5 p was evaluated within the 
p = 0 ~ 1.3 um range. To consider the impacts of both the shape and phase, we 
utilized the overlap integral n defined as 



E sim -E; ef ds 



§\E stm \ 2 ds§\ 



ds 



(3) 



Since rj was used primarily for very similar field profiles, the integration range was 
fixed to p = 0 ~ 7 um. 
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